Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε εδώ

1 Βασικές προϋποθέσεις

Εννοείται ότι πριν ξεκινήσουμε να κάνουμε το οτιδήποτε, έχουμε δημιουργήσει ένα νέο project στο R-Studio, το όνομα του οποίου - για λόγους τους οποίους θα καταλάβετε αργότερα - ΔΕΝ πρέπει να περιέχει:
1. ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Το ίδιο ισχύει ΚΑΙ για το file path του εν λόγω φακέλου (π.χ., όχι ‘E:/Εργαστήριο Οικολογίας/SDM’, αλλά ’E:/Ergastirio_Oikologias/SDM).

Ένα άλλο κρίσιμο σημείο είναι το εξής: θα χρειαστεί να εγκαταστήσετε τα πακέτα που πρόκειται να χρησιμοποιήσετε σε αυτό το tutorial1.

2 Περιβαλλοντικές Μεταβλητές

2.1 Εισαγωγή

Τα κλιματικά δεδομένα μας μπορούμε να τα κατεβάσουμε από το WorldClim, ενώ τα υψομετρικά δεδομένα από εδώ και τα δεδομένα σχετικά με την ξηρασία και την εξατμισιοδιαπνοή από εδώ. - Εάν φυσικά δουλεύουμε για μια συγκεκριμένη εργασία (πτυχιακή, μεταπτυχιακή, κτλ) και όχι στα πλαίσια ενός μαθήματος επιλογής. Προς το παρόν, θα χρησιμοποιήσουμε κάποιες εντολές για να κατεβάσουμε τα δεδομένα που θέλουμε μέσω της R.

Όπως μπορείτε να διαπιστώσετε, υπάρχουν πάρα πολλά κλιματικά μοντέλα (επί παραδείγματι, CCSM4, HadGEM-2, κτλ) και τέσσερα (4) διαφορετικά κλιματικά σενάρια που έχουν να κάνουν με τη συγκέντρωση των αερίων του θερμοκηπίου στην ατμόσφαιρα εντός του 21ου αιώνα { ενδεικτική βιβλιογραφία ή άλλη ενδεικτική βιβλιογραφία ή ακόμα μια ενδεικτική βιβλιογραφία ή τελευταία ενδεικτική βιβλιογραφία}2.

Συνεπώς, δεν μπορούμε να χρησιμοποιήσουμε μόνο ένα κλιματικό μοντέλο, πόσω μάλλον ένα κλιματικό σενάριο. Ευτυχώς όμως, δεν χρειάζεται και ούτε είναι επιστημονικά σωστό, να χρησιμοποιήσουμε όλα τα κλιματικά μοντέλα. Αναλόγως της περιοχής στην οποία εργαζόμαστε, κάποια κλιματικά μοντέλα αναπαριστούν καλύτερα τις περιβαλλοντικές συνθήκες σε σχέση με κάποια άλλα (McSweeney et al., 2015). Καθώς θα διαπιστώσετε (όταν με το καλό διαβάσετε την προαναφερθείσα εργασία), ακόμα και έτσι, είναι πολλά τα κλιματικά μοντέλα για την περιοχή μας. Όμως, μόνο 3 από αυτά έχουν δεδομένα και για τα 4 κλιματικά σενάρια (BCC, HadGEM, CCSM4). Εμείς χάριν συντομίας, θα εργαστούμε όμως μόνο με ένα (1) κλιματικό μοντέλο και ενα (1) κλιματικό σενάριο.

2.2 Προετοιμασία των μεταβλητών μας

2.2.1 Libraries

Ας φορτώσουμε τις απαραίτητες βιβλιοθήκες3.

library(easypackages)
libraries("sp", "rgdal", "raster", "dismo", "rgbif", "rgeos", "scales", "scrubr", 
    "mapr", "tidyverse", "rasterVis")

2.2.1.1 Download the raw data

Τώρα, ας δημιουργήσουμε τέσσερα αντικείμενα τα οποία θα περιέχουν τα αρχεία για τις παρούσες και τις μελλοντικές κλιματικές συνθήκες, τα υψομετρικά δεδομένα και την χώρα που μας ενδιαφέρει.

Πρώτα θα κατεβάσουμε τις γεωγραφικές μεταβλητές για την χώρα που μας ενδιαφέρει (οι ΗΠΑ στην προκειμένη περίπτωση). Θα πρέπει να εισάγουμε (στην εντολή country στον κώδικα που παρατίθεται πιο κάτω) τον 3-ψήφιο ISO κωδικό για τις ΗΠΑ που τον βρίσκουμε από αυτόν τον ιστότοπο ή με την εντολή ISO_countries <- getData("ISO3") %>% as.data.frame4. Εάν θέλουμε να περιορίσουμε την περιοχή μελέτης σε μια μικρότερη περιοχή, τότε θα πρέπει να χρησιμοποιήσουμε ένα άλλο λογισμικό, το Qgis, ώστε να “κόψουμε” την περιοχή που μας ενδιαφέρει. Προς το παρόν, δεν θα χρειαστεί να το κάνουμε αυτό, αλλά θα περιοριστούμε να κατεβάσουμε τα πλήρη δεδομένα για τις ΗΠΑ:

USA <- getData("GADM", country = "USA", level = 1)

Το αρχείο που φτιάξαμε περιέχει κάποια (αρκετά) δεδομένα. Ας τα δούμε και ας αναπαραστήσουμε γραφικά κάποια από αυτά:

datatable(USA@data)
Hawaii <- subset(USA, NAME_1 == "Hawaii")
SW_US <- subset(USA, NAME_1 == "California" | NAME_1 == "Nevada" | NAME_1 == 
    "Utah" | NAME_1 == "Arizona")

plot(USA)

plot(Hawaii)

plot(SW_US)

Κάντε το ίδιο για όποια πολιτεία της Αμερικής θέλετε. Στη συνέχεια, φτιάξτε ένα αντικείμενο το οποίο να περιέχει την Montana, το Idaho και το Wyoming και τέλος ένα αντικείμενο το οποίο να μην περιέχει την Αλάσκα και την Hawaii και ονομάστε το USA_ter.

Μετά θα κατεβάσουμε τις παρούσες και τις μελλοντικές κλιματικές συνθήκες5, όπως αυτές δίνονται από μια παγκόσμια διαδικτυακή βάση δεδομένων, το WorldClim. Οι κλιματικές αυτές μεταβλητές είναι 19 και είναι οι ακόλουθες:

Κλιματικές μεταβλητές

Var. Description
BIO1 Annual Mean Temperature
BIO2 Mean Diurnal Range (Mean of monthly (max temp – min temp))
BIO3 Isothermality (BIO2/BIO7) (* 100)
BIO4 Temperature Seasonality (standard deviation *100)
BIO5 Max Temperature of Warmest Month
BIO6 Min Temperature of Coldest Month
BIO7 Temperature Annual Range (BIO5-BIO6)
BIO8 Mean Temperature of Wettest Quarter
BIO9 Mean Temperature of Driest Quarter
BIO10 Mean Temperature of Warmest Quarter
BIO11 Mean Temperature of Coldest Quarter
BIO12 Annual Precipitation
BIO13 Precipitation of Wettest Month
BIO14 Precipitation of Driest Month
BIO15 Precipitation Seasonality (Coefficient of Variation)
BIO16 Precipitation of Wettest Quarter
BIO17 Precipitation of Driest Quarter
BIO18 Precipitation of Warmest Quarter
BIO19 Precipitation of Coldest Quarter

Όσον αφορά τις μελλοντικές κλιματικές συνθήκες, πρέπει να προσδιορίσουμε ποιό κλιματικό μοντέλο μας ενδιαφέρει (υπάρχουν πολλές επιλογές όπως προαναφέραμε, αλλά μας ενδιαφέρουν συγκεκριμένα κλιματικά μοντέλα), ποιό κλιματικό σενάριο (υπάρχουν 4 διαφορετικά: 26, 45, 60, 85), καθώς και ποιά χρονολογία μας ενδιαφέρει (50 ή 70 - αναφέρεται στα έτη 2050 και 2070, αντίστοιχα). Προς το παρόν εμείς θα δουλέψουμε με το μοντέλο CCMS4 (αντιστοιχεί στο μοντέλο CC), το κλιματικό σενάριο 85 και το έτος 2070. Τις μελλοντικές κλιματικές συνθήκες θα τις χρησιμοποιήσουμε προκειμένου να δείτε πώς θα προβλέψουμε την κατανομή των ειδών στον χώρο σε άλλο χρόνο:

current_climate <- getData("worldclim", var = "bio", res = 2.5)

current_climate
## class       : RasterStack 
## dimensions  : 3600, 8640, 31104000, 19  (nrow, ncol, ncell, nlayers)
## resolution  : 0.04166667, 0.04166667  (x, y)
## extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names       :  bio1,  bio2,  bio3,  bio4,  bio5,  bio6,  bio7,  bio8,  bio9, bio10, bio11, bio12, bio13, bio14, bio15, ... 
## min values  :  -278,     9,     8,    64,   -86,  -559,    53,  -278,  -501,  -127,  -506,     0,     0,     0,     0, ... 
## max values  :   319,   213,    96, 22704,   489,   258,   725,   376,   365,   382,   289, 10577,  2437,   697,   265, ...
future_climate <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, model = "CC", 
    year = 70)

future_climate
## class       : RasterStack 
## dimensions  : 3600, 8640, 31104000, 19  (nrow, ncol, ncell, nlayers)
## resolution  : 0.04166667, 0.04166667  (x, y)
## extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
## names       : cc85bi701, cc85bi702, cc85bi703, cc85bi704, cc85bi705, cc85bi706, cc85bi707, cc85bi708, cc85bi709, cc85bi7010, cc85bi7011, cc85bi7012, cc85bi7013, cc85bi7014, cc85bi7015, ... 
## min values  :      -217,         9,         7,        54,       -58,      -491,        53,      -284,      -436,        -94,       -443,          0,          0,          0,          0, ... 
## max values  :       351,       219,        96,     21959,       536,       281,       707,       426,       411,        430,        325,      13161,       3017,        682,        223, ...

Ας αλλάξουμε τα ονόματα των μεταβλητών μας ως ακολούθως:

names(current_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

names(future_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

Καθώς τα κλιματικά μας δεδομένα όπως διαπιστώσατε έχουν εξαιρετικά μεγάλο εύρος τιμών, θα χρειαστεί να τα μετατρέψουμε σε oC όσον αφορά την θερμοκρασία (bio_1). Εδώ μπορείτε να διαπιστώσετε το γιατί6. Αυτό γίνεται με την ακόλουθη εντολή:

gain(current_climate$bio_1) = 0.1

gain(future_climate$bio_1) = 0.1

Στη συνέχεια, θα κατεβάσουμε τα υψομετρικά δεδομένα για την χώρα που μας ενδιαφέρει:

altitude <- getData("alt", country = "USA", mask = TRUE)

altitude
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA1_msk_alt.grd`
## class       : RasterLayer 
## dimensions  : 3012, 6948, 20927376  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -124.8, -66.9, 24.4, 49.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA1_msk_alt.grd 
## names       : USA1_msk_alt 
## values      : -171, 4275  (min, max)
## 
## 
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA2_msk_alt.grd`
## class       : RasterLayer 
## dimensions  : 2448, 5916, 14482368  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -179.2, -129.9, 51.1, 71.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA2_msk_alt.grd 
## names       : USA2_msk_alt 
## values      : -22, 6098  (min, max)
## 
## 
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA3_msk_alt.grd`
## class       : RasterLayer 
## dimensions  : 228, 900, 205200  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 172.4, 179.9, 51.2, 53.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA3_msk_alt.grd 
## names       : USA3_msk_alt 
## values      : -1, 1099  (min, max)
## 
## 
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA4_msk_alt.grd`
## class       : RasterLayer 
## dimensions  : 420, 672, 282240  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -160.3, -154.7, 18.8, 22.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA4_msk_alt.grd 
## names       : USA4_msk_alt 
## values      : -3, 4120  (min, max)
## We are only interested in the first of the four elements of the
## altitudinal data

altitude <- altitude[[1]]

2.2.1.2 Crop and mask

Στη συνέχεια, θα ‘κόψουμε’ τις περιβαλλοντικές μας μεταβλητές στο σχήμα των ΗΠΑ και ακολούθως θα δημιουργήσουμε ένα νέο αρχείο που θα περιέχει τις μεταβλητές αυτές και την περιοχή μελέτης μας:

cropped_current <- crop(current_climate, USA_ter, snap = "in") %>% mask(., USA_ter)

cropped_future <- crop(future_climate, USA_ter, snap = "in") %>% mask(., USA_ter)

Τώρα, θα ‘κόψουμε’ τις κλιματικές μεταβλητές για τις ΗΠΑ σε επιμέρους αρχεία:

SW_crop_current <- crop(cropped_current, SW_US, snap = "in") %>% mask(., SW_US)

SW_crop_future <- crop(cropped_future, SW_US, snap = "in") %>% mask(., SW_US)

Κάντε το ίδιο για την California και το Oregon.

Μετά, μπορούμε να αναπαραστήσουμε γραφικά τα αποτελέσματα:

gplot(cropped_current[[1]]) + geom_raster(aes(fill = value)) + scale_fill_gradientn(colours = c("brown", 
    "red", "yellow", "darkgreen", "green")) + coord_equal() + xlab("Longitude") + 
    ylab("Latitude") + ggtitle("Temperature")

Κάντε το ίδιο για το αρχείο που φτιάξατε προηγουμένως.

Στη συνέχεια, θα υπολογίσουμε κάποια βασικά στατιστικά στοιχεία για κάθε αρχείο το οποίο δημιουργήσαμε:

cellStats(cropped_current, "quantile")
##      bio_1 bio_2 bio_3 bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10 bio_11
## 0%    -7.1    59    22  1794    60  -245   120  -135  -157      5   -157
## 25%    6.6   123    31  7431   278  -130   350   105   -44    184    -56
## 50%    9.8   134    37  8385   305   -85   382   175    36    213    -13
## 75%   14.6   155    41  9626   329   -28   419   211   161    248     46
## 100%  25.2   213    65 13060   456   178   509   332   329    360    212
##      bio_12 bio_13 bio_14 bio_15 bio_16 bio_17 bio_18 bio_19
## 0%       46      7      0      5     19      0      1     11
## 25%     392     66     11     20    170     39    121     51
## 50%     692     96     21     36    262     77    212    112
## 75%    1067    121     53     56    331    187    296    233
## 100%   3345    556    154    104   1512    480    660   1455
cellStats(cropped_current, "mean")
##      bio_1      bio_2      bio_3      bio_4      bio_5      bio_6 
##   10.59007  138.48284   36.36840 8517.96950  301.83563  -79.58974 
##      bio_7      bio_8      bio_9     bio_10     bio_11     bio_12 
##  381.42537  153.05698   55.71655  213.60417   -5.79786  747.70502 
##     bio_13     bio_14     bio_15     bio_16     bio_17     bio_18 
##   98.39876   32.30871   38.61828  267.33779  113.08449  209.89919 
##     bio_19 
##  161.71277
cellStats(cropped_current, "sd")
##       bio_1       bio_2       bio_3       bio_4       bio_5       bio_6 
##    5.268230   21.448783    6.533989 1729.940968   38.814827   69.275880 
##       bio_7       bio_8       bio_9      bio_10      bio_11      bio_12 
##   54.956739   83.319657  113.217157   44.299732   67.841672  415.811072 
##      bio_13      bio_14      bio_15      bio_16      bio_17      bio_18 
##   51.187839   26.819703   20.246943  144.910739   88.540570  112.976747 
##      bio_19 
##  151.904632

Κάντε το ίδιο για το SW_crop και για το αρχείο που φτιάξατε για την California και το Oregon.

Τέλος, μπορούμε να δημιουργήσουμε υποσύνολα των περιοχών που μας ενδιαφέρουν, ώστε να ανταποκρίνονται σε συγκεκριμένες κλιματικές (ή άλλες συνθήκες):

hot_US <- cropped_current[[1]] < 25 & cropped_current[[1]] > 14.6

plot(hot_US)

Δοκιμάστε το και για άλλες περιοχές και μεταβλητές.

2.2.1.3 Σύγκριση

Τώρα, θα πρέπει να ελέγξουμε εάν τα υψομετρικά δεδομένα έχουν την ίδια ανάλυση (resolution), μέγεθος (nrow, ncol) και έκταση (bbox), με τα περιβαλλοντικά μας δεδομένα:

altitude
## class       : RasterLayer 
## dimensions  : 3012, 6948, 20927376  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -124.8, -66.9, 24.4, 49.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA1_msk_alt.grd 
## names       : USA1_msk_alt 
## values      : -171, 4275  (min, max)
cropped_current
## class       : RasterBrick 
## dimensions  : 596, 1387, 826652, 19  (nrow, ncol, ncell, nlayers)
## resolution  : 0.04166667, 0.04166667  (x, y)
## extent      : -124.75, -66.95833, 24.54167, 49.375  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       :   bio_1,   bio_2,   bio_3,   bio_4,   bio_5,   bio_6,   bio_7,   bio_8,   bio_9,  bio_10,  bio_11,  bio_12,  bio_13,  bio_14,  bio_15, ... 
## min values  :    -7.1,    59.0,    22.0,  1794.0,    60.0,  -245.0,   120.0,  -135.0,  -157.0,     5.0,  -157.0,    46.0,     7.0,     0.0,     5.0, ... 
## max values  :    25.2,   213.0,    65.0, 13060.0,   456.0,   178.0,   509.0,   332.0,   329.0,   360.0,   212.0,  3345.0,   556.0,   154.0,   104.0, ...

2.2.1.4 Resolution change

Στην περίπτωση κατά την οποία τα δύο σετ δεδομένων έχουν διαφορετική ανάλυση, έκταση και μέγεθος, θα χρειαστεί να αλλάξουμε την ανάλυση των υψομετρικών δεδομένων7.

agg_altitude <- aggregate(altitude, 
                 5, ## Why choose 5 and not any other number?
                 fun = mean)

2.2.1.5 Pixel size

Μετά, θα αλλάξουμε το μέγεθος των pixel των υψομετρικών δεδομένων για να συμφωνεί με αυτό των περιβαλλοντικών δεδομένων:

elev_US_resample <- resample(agg_altitude, cropped_current, method = "bilinear")

2.2.1.6 Final Stack

Τέλος, θα τα ενώσουμε σε ένα αρχείο, τόσο για το παρόν, όσο και για το μέλλον:

var_US <- readAll(stack(cropped_current, elev_US_resample))

var_US_70 <- readAll(stack(cropped_future, elev_US_resample))


elev_SW_US <- crop(elev_US_resample, SW_US, snap = "in") %>% mask(., SW_US)

var_SW_US <- readAll(stack(SW_crop_current, elev_SW_US))

var_SW_70 <- readAll(stack(SW_crop_future, elev_SW_US))

Καλό είναι να αποθηκεύσετε τα αρχεία αυτά υπό την μορφή Rds και να τα αφαιρέσετε από την μνήμη της R8. Αυτό το κάνουμε, προκειμένου να μπορούμε εύκολα να το φορτώνουμε στην R και να μην υπάρχουν τα δεδομένα αυτά συνεχώς στην μνήμη της session που έχουμε ανοικτή, διότι αυξάνουν αρκετά τον χρόνο απόκρισης της (καθυστερεί πολύ η R). Μπορείτε ανά πάσα στιγμή να τα φορτώσετε στην R, με την ακόλουθη εντολή:

## Creating a folder to save our data
dir.create("./Data/Variables", recursive = TRUE, showWarnings = FALSE)

## Save our data in rds format
saveRDS(var_US, "./Data/Variables/Environmental variables and Elevation (30s) USA.rds")
saveRDS(var_US_70, "./Data/Variables/Environmental variables and Elevation (30s) USA 2070.rds")

## Load our data from rds format
var_US <- readRDS("./Data/Variables/Environmental variables and Elevation (30s) USA.rds")
var_US_70 <- readRDS("./Data/Variables/Environmental variables and Elevation (30s) USA 2070.rds")

2.2.1.7 Plot

Ας αναπαραστήσουμε γραφικά το αποτέλεσμα, το οποίο θα φανεί στην ακόλουθη εικόνα:

plot(var_SW_US[[c(1, 12, 20)]], col = rainbow(100, start = 0, end = 0.8))

Τρέξτε την εντολή cellStats για τα 3 αρχεία τα οποία έχετε φτιάξει (var_US, var_SW_US και var_CO_US). Βλέπετε διαφορές;

Ας αναπαραστήσουμε γραφικά - αλλά με τέτοιον τρόπο ώστε να μας δίνει κάποιες στατιστικές πληροφορίες - τις περιβαλλοντικές μας μεταβλητές:

histogram(var_SW_US[[20]])

levelplot(var_SW_US[[20]])

Κάντε το ίδιο και για άλλες μεταβλητές

Ας δούμε τα δεδομένα μας τρισδιάστατα:

plot3D(var_SW_US[[20]])

Κάντε το ίδιο και για άλλες μεταβλητές

2.3 Συσχέτιση των μεταβλητών μας

Σε αυτήν την ενότητα θα διαπιστώσουμε εάν οι περιβαλλοντικές μας μεταβλητές έχουν υψηλή συσχέτιση μεταξύ τους στην περιοχή μελέτης.

2.3.1 Φόρτωση των απαραίτητων βιβλιοθηκών

library(usdm)

2.3.2 Μετατροπή των raster αρχείων σε data-frames

Στο σημείο αυτό θα μετατρέψουμε τα raster αρχεία μας σε data-frames, προκειμένου να μπορέσουμε να διαπιστώσουμε ποιές μεταβλητές εμφανίζουν υψηλή συσχέτιση:

var_US_df <- na.omit(as.data.frame(var_US))

var_SW_US_df <- na.omit(as.data.frame(var_SW_US))

Γιατί το κάνουμε αυτό μόνο για τις παρούσες κλιματικές συνθήκες;

2.3.3 Έλεγχος συγγραμικότητας και Variance Inflation Factor (VIF)

Όπως θα γνωρίζετε, ΔΕΝ μπορούμε να χρησιμοποιήσουμε σε μια ανάλυση δύο μεταβλητές, οι οποίες εμφανίζουν συντελεστή συσχέτισης \(> 0.7\)9, καθότι στην προκειμένη περίπτωση οι δύο αυτές μεταβλητές είναι συγγραμικές και ως εκ τούτου εάν τις συμπεριλάβουμε και τις δύο στην ανάλυση μας, τότε δεν θα γνωρίζουμε σε ποιά εκ των δύο θα οφείλεται το αποτέλεσμα της ανάλυσης μας.

Προκειμένου να ξεπεράσουμε αυτόν τον σκόπελο, θα τρέξουμε την ακόλουθη εντολή10:

vifcor(var_SW_US_df, th = 0.7)
## 13 variables from the 20 input variables have collinearity problem: 
##  
## bio_16 bio_6 bio_19 bio_12 bio_10 bio_17 bio_11 bio_7 bio_3 bio_5 bio_1 bio_14 bio_15 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( bio_18 ~ bio_13 ):  0.02776225 
## max correlation ( bio_13 ~ bio_4 ):  -0.6901746 
## 
## ---------- VIFs of the remained variables -------- 
##      Variables      VIF
## 1        bio_2 1.281558
## 2        bio_4 2.403638
## 3        bio_8 2.160976
## 4        bio_9 1.807053
## 5       bio_13 2.232166
## 6       bio_18 2.482672
## 7 USA1_msk_alt 4.043632

Τώρα δοκιμάστε το και εσείς, αλλάζοντας το όριο σε 0.8 και 0.9. Τι παρατηρείτε; Ποιές μεταβλητές θα κρατούσατε στην ανάλυση;

2.3.4 Επιλογή των μεταβλητών μας

Πλεόν είμαστε σε θέση να επιλέξουμε τις μη συγγραμικές μεταβλητές, τις οποίες θα χρησιμοποιήσουμε στην ανάλυση μας.

## If these are our uncorrelated variables
var_SW_US_sub <- stack(subset(var_SW_US, c(2, 3, 8, 13:15, 20)))
var_SW_70_sub <- stack(subset(var_SW_70, c(2, 3, 8, 13:15, 20)))

## Let's scale our variables (mean = 0, sd = 1)
var_SW_US_sub <- stack(scale(var_SW_US_sub))
var_SW_70_sub <- stack(scale(var_SW_70_sub))

Κάντε το ίδιο για το σύνολο των ΗΠΑ.

3 Δεδομένα θέσης

3.1 Εισαγωγή

Στο σημείο αυτό, είμαστε πλέον στην ευχάριστη θέση να μπορούμε να χρησιμοποιήσουμε τα δεδομένα θέσης για τα είδη τα οποία μας ενδιαφέρουν. Υπάρχουν δύο περιπτώσεις, όσον αφορά την απόκτηση τέτοιων δεδομένων:
- να τα έχουμε αποκτήσει πρωτογενώς (δηλαδή να έχουμε εργαστεί στο πεδίο και να έχουμε εντοπίσεις τις θέσεις εμφάνισης του εκάστοτε είδους και των υποπληθυσμών του)
- να τα έχουμε αποκτήσει δευτερογενώς (δηλαδή είτε να έχουμε χρησιμοποιήσει τις διαδικτυακές βάσεις ψηφιοποιημένων - και ενίοτε γεωαναφερμένων - δεδομένων GBIF και BIEN, είτε να έχουμε ψηφιοποιήσει χάρτες εμφάνισης θέσης11)

Το ιδανικό βεβαίως είναι να έχουμε αποκτήσει πρωτογενώς τα στοιχεία αυτά, καθότι στην περίπτωση αυτή, μπορούμε να είμαστε σίγουροι για τις γεωγραφικές συντεταγμένες των υποπληθυσμών του προς μελέτη είδους. Στην πλειονότητα των περιπτώσεων όμως, δεν έχουμε πρόσβαση σε τέτοια δεδομένα, οπότε αναγκαστικά καταφεύγουμε στην δεύτερη λύση και δη στις διαδικτυακές βάσεις δεδομένων. Ευτυχώς για εμάς, υπάρχουν κάποια πακέτα12 στην R, όπως το rgbif και το rbien, τα οποία αυτοματοποιούν την διαδικασία μεταφόρτωσης δεδομένων από τις διαδικτυακές αυτές βάσεις13.

Εμείς θα εργαστούμε με το πακέτο rgbif προς το παρόν.

3.2 Φόρτωση των απαραίτητων βιβλιοθηκών

library(easypackages)
libraries("viridis", "ggmap", "biogeo", "spThin", "mapr")

3.3 Ονοματολογικός έλεγχος

Σήμερα θα ασχοληθούμε με ένα είδος το οποίο είναι ενδημικό της δυτικής ακτής των ΗΠΑ, το Yucca brevifolia.

Το είδος Yucca brevifolia στην πολιτεία της Nevada

Το είδος Yucca brevifolia στην πολιτεία της Nevada

Σύμφωνα με το Υπουργείο Αγροτικής Ανάπτυξης των ΗΠΑ, το είδος αυτό είναι ενδημικό των πολιτειών της California, της Nevada, της Utah και της Arizona. Ανήκει στην οικογένεια Agavaceae και το κοινό του όνομα είναι Joshua tree, εξαιτίας της μορφής του: οι Μορμόνοι άποικοι οι οποίοι αντίκρυσαν για πρώτη φορά αυτό το δένδρο στην έρημο Mojave τον 18ο αιώνα θέωρησαν ότι τους θύμιζε τον ικέτη Ιησού που σηκώνει τα χέρια του προς τον ουρανό καθώς προσεύχεται.

Ο Άγιος Ιησούς ο Δίκαιος μαζί με τον Μωϋσή στην Καναάν (Πίνακας του Giovanni Lafranco, από το J. Paul Getty Museum του Los Angeles)
Ένα άλλο κοινό όνομα του είδους αυτού είναι το izote de desierto (Ισπανικά: το μαχαίρι της ερήμου). Η Yucca brevifolia εντοπίζεται κυρίως στην έρημο Mojave σε υψόμετρο μεταξύ 400 έως και 1800 m, είναι ένα ταχυαυξές είδος (7.6 cm/έτος για τα πρώτα 10 χρόνια) το οποίο δεν σχηματίζει ετήσιους δακτύλιους. Το ρίζικο του σύστημα φθάνει σε βάθος 11 m. Τα ψηλότερα άτομα έχουν 15 m ύψος, ενώ η περίοδος ανθοφορίας του είναι μεταξύ Φεβρουαρίου και Απριλίου: όπως τα περισσότερα είδη τα οποία απαντώνται σε ερήμους, χρειάζονται τη σωστή ποσότητα βροχής το κατάλληλο χρονικό διάστημα, προκειμένου να ανθίσουν. Έχουν δε, βρεθεί άτομα ηλικίας 1000 ετών. Σύμφωνα με τους Shafer et al. (2001),η κλιματική αλλαγή αναμένεται να έχει αρνητική επίδραση στο είδος αυτό.

Πρώτα χρειάζεται να βεβαιωθούμε ότι το προς μελέτη είδος, υπάρχει όντως εντός της διαδικτυακής βάσης την οποία επιθυμούμε να χρησιμοποιήσουμε, προκειμένου να εξάγουμε δεδομένα. Η εντολή name_suggest() ερευνά όλα τα είδη (taxa καλύτερα) τα οποία μοιάζουν με το είδος το οποίο ψάχνουμε και κρατά μόνο το όνομα των taxa το οποίο ταιριάζει ακριβώς με αυτό το οποίο μας ενδιαφέρει14.

spp_int <- name_suggest(q = "Yucca brevifolia", rank = "species", limit = 10000)
(spp_int <- spp_int[grepl("^Yucca brevifolia$", spp_int$canonicalName), ])
key canonicalName rank
7520449 Yucca brevifolia SPECIES
7636266 Yucca brevifolia SPECIES
2775592 Yucca brevifolia SPECIES
7390928 Yucca brevifolia SPECIES

Στο βήμα αυτό αναγνωρίσαμε τον μοναδικό αριθμό αναφοράς (key) του προς μελέτη είδους και μόνο χρησιμοποιώντας τον αριθμό αυτό μπορούμε να είμαστε βέβαιοι ότι τα δεδομένα θέσης τα οποία θα λάβουμε από την διαδικτυακή αυτή βάση δεδομένων θα αναφέρονται στο είδος το οποίο μας ενδιαφέρει.

3.4 Μεταφόρτωση των γεωαναφερμένων δεδομένων θέσης

sp_occ <- occ_search(taxonKey = spp_int$key, country = "US", fields = c("name", 
    "key", "country", "decimalLatitude", "decimalLongitude"), hasCoordinate = T, 
    limit = 1000, return = "data")

sp_occ
## $`7520449`
## [1] "no data found, try a different search"
## 
## $`7636266`
## [1] "no data found, try a different search"
## 
## $`2775592`
## # A tibble: 1,000 x 5
##    name                    key decimalLongitude decimalLatitude country   
##    <chr>                 <int>            <dbl>           <dbl> <chr>     
##  1 Yucca brevifolia 1806362560             -116            33.9 United St~
##  2 Yucca brevifolia 1807289573             -116            34.1 United St~
##  3 Yucca brevifolia 1807298279             -114            36.6 United St~
##  4 Yucca brevifolia 1807303045             -116            34.1 United St~
##  5 Yucca brevifolia 1806341722             -115            36.2 United St~
##  6 Yucca brevifolia 1805394434             -116            34.1 United St~
##  7 Yucca brevifolia 1831095189             -116            34.0 United St~
##  8 Yucca brevifolia 1805399185             -116            34.0 United St~
##  9 Yucca brevifolia 1805417949             -114            36.7 United St~
## 10 Yucca brevifolia 1831185811             -114            36.5 United St~
## # ... with 990 more rows
## 
## $`7390928`
## [1] "no data found, try a different search"
## 
## attr(,"args")
## attr(,"args")$hasCoordinate
## [1] TRUE
## 
## attr(,"args")$limit
## [1] 1000
## 
## attr(,"args")$offset
## [1] 0
## 
## attr(,"args")$taxonKey
## [1] 7520449 7636266 2775592 7390928
## 
## attr(,"args")$country
## [1] "US"
## 
## attr(,"args")$fields
## [1] "name"             "key"              "country"         
## [4] "decimalLatitude"  "decimalLongitude"

Προκειμένου να μην προκύψουν προβλήματα με το μονοπάτι (path) του φακέλου και των δεδομένων μας, καλό είναι να αφαιρέσουμε τα κενά εντός του ονόματος του είδους μας και να διαλέξουμε το ID15 του taxon μας με τις περισσότερες εγγραφές:

sp_occ <- sp_occ[["2775592"]]
sp_occ$name <- sub(" ", ".", sp_occ$name)
(spp_to_model <- unique(sp_occ$name))
## [1] "Yucca.brevifolia"

Τώρα έχει έρθει η στιγμή να δούμε πόσα γεωαναφερμένα δεδομένα θέσης έχουμε στην διάθεση μας:

sort(table(sp_occ$name), decreasing = T)
## Yucca.brevifolia 
##             1000

Ας οπτικοποιήσουμε τα αποτελέσματα της ανάλυσης μας έως αυτό το σημείο:

## Setting the extent (data taken from our RasterStack)
myloc <- c(-125, -67, 24.6, 49.3)

mymap <- get_map(location = myloc,
                 source = 'stamen',
                 maptype = 'terrain',
                 crop = F)

## Vizualise our map
ggmap(mymap) + 
  geom_point(aes(x = decimalLongitude, y = decimalLatitude), 
             data = sp_occ, 
             alpha = .5, 
             color = "darkred", 
             size = 3) + 
  xlab('Longitude') +
  ylab('Latitude')

Θα πρέπει να τονίσουμε το εξής: Εφ’όσον δεν έχουμε συλλέξει πρωτογενώς τα δεδομένα μας, θα χρειαστεί να ελέγξουμε εάν οι θέσεις εμφάνισης που μεταφορτώσαμε από τις διαδικτυακές βάσεις δεδομένων, συμφωνούν με τα γνωστά όρια εξάπλωσης του είδους που μας ενδιαφέρει. Εάν όντως υπάρχουν κάποια σημεία τα οποία βρίσκονται εκτός της περιοχής εξάπλωσης, θα χρειαστεί να τα αφαιρέσουμε από την ανάλυση. Θα χρησιμοποιήσουμε το πακέτο scrubr για να κάνουμε αυτή την ανάλυση16.

## Let's create a new object, containing only the correct coords
sp_occ_ck <- dframe(sp_occ) %>% coord_impossible() %>% coord_incomplete() %>% 
    coord_unlikely()

## Is there any difference between the two sets of coords?
NROW(sp_occ)
## [1] 1000
NROW(sp_occ_ck)
## [1] 1000

Ας φτιάξουμε έναν φάκελο ώστε να αποθηκεύσουμε τα δεδομένα μας και ας τα σώσουμε, προκειμένου να μπορέσουμε να τα φορτώσουμε κάποια άλλη στιγμή:

dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Yucca brevifolia", recursive = TRUE, showWarnings = FALSE)

saveRDS(sp_occ_ck, "./Species/Yucca brevifolia.rds")

4 Δεδομένα θέσης και περιβαλλοντικές μεταβλητές

4.1 Εισαγωγή

Πλέον μπορούμε να συνδέσουμε τις γεωαναφερμένες θέσεις εμφάνισης του είδους που μας ενδιαφέρει με τα raster αρχεία τα οποία είχαμε δημιουργήσει προηγουμένως. Στο σημείο αυτό, μπορούμε να εξάγουμε τα περιβαλλοντικά δεδομένα για κάθε μια εκ των θέσεων εμφάνισης του είδους που μας ενδιαφέρει με μια μόνο εντολή, αναδεικνύοντας κατ’αυτόν τον τρόπο το πλεονέκτημα του να έχουμε ενώσει προηγουμένως τα raster αρχεία μας σε ένα RasterStack.

coordinates(sp_occ_ck) <- c("longitude", "latitude")

pts.clim <- raster::extract(var_SW_US, sp_occ_ck, method = "bilinear")

sp_occ_clim <- data.frame(cbind(coordinates(sp_occ_ck), pts.clim, sp_occ_ck@data))

Ακολούθως, θέλουμε να αφαιρέσουμε τις θέσεις εμφάνισης οι οποίες πέφτουν εκτός των περιβαλλοντικών μας μεταβλητών (μπορεί να πέφτουν στην θάλασσα).

sp_occ_clim <- sp_occ_clim %>% na.omit()
coordinates(sp_occ_clim) <- c("longitude", "latitude")

Και τώρα, ήρθε η στιγμή να οπτικοποιήσουμε τα αποτελέσματα μας και δη, την παρουσία του είδους μας στην περιοχή μελέτης:

map.ext <- extent(var_SW_US)

plot(var_SW_US[[20]], col = terrain.colors(100), ext = map.ext)

plot(sp_occ_clim, pch = 16, cex = 0.5, add = T)

plot(SW_US, add = T)

Στο σημείο αυτό, θα ήθελα να σας τονίσω ότι ορισμένες φορές, πέρα από τους ελέγχους που διενεργήσαμε έως τώρα, ίσως χρειαστεί να προχωρήσουμε σε μια αραίωση (thinning) των θέσεων εμφάνισης. Αυτό μπορεί να είναι απαραίτητο όταν έχουμε πάρα πολλά σημεία τα οποία είναι αρκετά κοντά το ένα στο άλλο17. Σε κάθε περίπτωση, το επιστημονικά ορθό είναι να υπάρχει μια (1) θέση εμφάνισης ανά km2. Ένα πακέτο που μπορεί να μας βοηθήσει είναι το spThin18.

## Min. distance: 1000 m Number of reps: 1000
spoccs <- correct_occurrence_points_sp %>% as.data.frame() %>% rename(x = coords.x1, 
    y = coords.x2, Species = name) %>% addmainfields(species = "Species")

thinned <- thin(loc.data = spoccs, lat.col = "y", long.col = "x", spec.col = "Species", 
    thin.par = 1, reps = 1000, verbose = T, out.dir = "./Thinned", locs.thinned.list.return = TRUE, 
    write.files = TRUE, out.base = "Yucca brevifolia thinned")

Βιβλιογραφία

McSweeney, C.F., Jones, R.G., Lee, R.W., & Rowell, D.P. (2015) Selecting CMIP5 GCMs for downscaling over multiple regions. Clim. Dyn., 44, 3237–3260.

Shafer, S.L., Bartlein, P.J., & Thompson, R.S. (2001) Potential changes in the distributions of western north america tree and shrub taxa under future climate scenarios. Ecosystems, 4, 200–215.

  1. Υπάρχουν δύο τρόποι για να το κάνετε αυτό γρήγορα:
    1. library(easypackages) packages('package_name_1', 'package_name_2')
    2. install.packages(c('package_name_1', 'package_name_2'), dependecies = TRUE)↩

  2. Καλό είναι να τις διαβάσετε↩

  3. Εάν δεν τις έχετε εγκαταστήσει, κάντε το τώρα: install.packages(insert_package_name_here) ή library(easypackages), packages('raster', 'sp', 'rgdal') - εάν θέλετε να εγκαταστήσετε τα τρια αυτά πακέτα↩

  4. Ή με αυτήν την εντολή εάν θέλουμε μόνο τις ΗΠΑ: ISO_countries <- getData("ISO3") %>% as.data.frame %>% filter(NAME=="USA")↩

  5. Επειδή θα πάρει αρκετό χρόνο με τη σύνδεση στο εργαστήριο, να τις έχετε κατεβάσει. Διαφορετικά, ζητήστε μου να σας τις δώσω εγώ↩

  6. Γιατί;↩

  7. Γιατί όμως να μην αλλάξουμε την ανάλυση των περιβαλλοντικών δεδομένων;↩

  8. Πώς τα αφαιρούμε από την μνήμη της R;↩

  9. Αυτό είναι το λεγόμενο hard boundary - άλλοι ερευνητές εφαρμόζουν το όριο \(> 0.8\) ή ακόμα και \(> 0.9\), αλλά είναι επιστημονικά ασφαλέστερο να εφαρμόζει κανείς το αυστηρότερο όριο↩

  10. Εδώ εφαρμόζουμε το αυστηρότερο όριο - μπορείτε να το αλλάξετε, εάν επιθυμείτε↩

  11. Ένα τυπικό παράδειγμα είναι ο Άτλαντας της Χλωρίδας του Αιγαίου ή ο Άτλαντας της Χλωρίδας της Ευρώπης, όμως στην περίπτωση αυτή υπεισέρχονται άλλα προβλήματα τα οποία σχετίζονται με την ανάλυση των εν λόγω χαρτών (π.χ., ο Άτλαντας της Χλωρίδας της Ευρώπης είναι σε ανάλυση \(50 \times 50 km^2\))↩

  12. Δεν είναι τα μοναδικά όμως: υπάρχουν πολλά άλλα, όπως το dismo και το spocc↩

  13. Κάθε ένα από αυτά τα πακέτα, έχει τουλάχιστον μια σύντομη περιγραφή (αγγλιστί: vignette), την οποία είναι ΑΠΑΡΑΙΤΗΤΟ να διαβάσετε - ειδικά αυτές του rgbif↩

  14. Με την εντολή ?rgbif::name_suggest() μπορείτε να δείτε τα arguments που δέχεται η εντολή αυτή και τι μπορείτε να αλλάξετε - σας συνιστώ εντόνως να το κάνετε, αλλά και για κάθε εντολή που χρησιμοποιείτε↩

  15. Εάν υπάρχουν πολλά βέβαια↩

  16. Δείτε το manual του πακέτου scrubr ή τρέξτε τις εντολές μια μια για να δείτε τι κάνουν στον κώδικα που έχετε στα χέρια σας↩

  17. Σε απόσταση \(<1000m\) όταν εργαζόμαστε με χάρτες υποβάθρου κλίμακας 1 x 1 km2 (ή ανάλυσης 30 arc-seconds)↩

  18. Θα σας βοηθήσει το vignette του εν λόγου πακέτου↩